www.gusucode.com > 循环自相关函数工具箱源码程序 > matlab代做 修改 程序循环自相关函数工具箱/cyclostationary_toolbox/cyclic_cross_periodogram.m
function S=cyclic_cross_periodogram(x,y,N,a,g,L) % % CYCLIC_CROSS_PERIODOGRAM % % calculate the spectral correlation density function % of signals x and y using the cyclic cross periodogram % algorithm % % Input parameters are: % x,y signals % N length of time window used for estimating frequency % segments % a window used for smoothing segments % g window for smoothing correlation % L decimation factor % % USAGE S=cyclic_cross_periodogram(x,y,N,a,g,L) % % e.g s=cyclic_cross_periodogram(s1,s1,64,'hamming','hamming',1) if ~exist('L') L=1; end lx=length(x); if (length(y)~=lx) error('Time series x and y must be same length') end n=0:floor((lx-N)/L); ln=length(n); a=feval(a,N)'; g=feval(g,ln)'; g=g/sum(g); a=a/sum(a); X=zeros(2*N+1,ln); Y=zeros(2*N+1,ln); Ts=1/N; for f=-N:N xf=x.*exp(-j*2*pi*f*(0:lx-1)*Ts); yf=y.*exp(-j*2*pi*f*(0:lx-1)*Ts); for i=1:ln n_r=n(i)*L+(1:N); X(f+N+1,i)=sum(a.*xf(n_r)); Y(f+N+1,i)=conj(sum(a.*yf(n_r))); end end S=zeros(N+1,floor(N/2)+1); for alpha=-floor(N/4):floor(N/4) for f=-floor(N/2):floor(N/2) f1=f+alpha; f2=f-alpha; if ((abs(f1)<N/2)&(abs(f2)<N/2)) S(f+floor(N/2)+1,floor(N/4)+alpha+1)=sum(g.*X(f1+N+1,:).*Y(f2+N+1,:)); end end end